
void psi_angle(){


	const double Pi = TMath::Pi();
	const double N_particle = 10000;
	
	double Psi[4] = {0,0,0,0};

	TRandom *myRn = new TRandom();
	for (int ih = 1; ih<4; ih++){
			Psi[ih] = (myRn->Uniform(-1*Pi / ih , Pi / ih));
							}
	double v2_input = 0.5;

	TF1 *dNdphi = new TF1("dNdphi","[0] *(   [1]+  2*[2]*cos(1*(x-[3])) +2*     [4]*cos(2*(x-[5])) +2* [6]*cos(3*(x-[7])))", 0,2* Pi);
	dNdphi->SetParameter(0, 100);    //coeffi
	dNdphi->SetParameter(1, 1);       // const term
	dNdphi->SetParameter(2, 0);      // directed flow magnitude  (v1)
	dNdphi->SetParameter(3, Psi[1]);     // directed flow direction  (psi 1)
	dNdphi->SetParameter(4, v2_input);   // elliptic flow magnitude (v2) 
	dNdphi->SetParameter(5, Psi[2]);     // elliptic flow direction (psi 2) 
	dNdphi->SetParameter(6, 0);      // tri flow magnitude (v3)
	dNdphi->SetParameter(7, Psi[3]);     // tri flow direction (psi 3)


	TCanvas *c1 = new TCanvas();
	TH1F *h1 = new TH1F("phi dist", "phi dist", 4400, 0, Pi);
	TH1F *h2 = new TH1F("","",100,0,Pi);
	TGraph *g1 = new TGraph();
	c1->Divide(1,2);

	double phi[N_particle];
	for (int i_particle=0; i_particle<N_particle; i_particle++){
		phi[i_particle] = dNdphi->GetRandom() ;
		h1->Fill(phi[i_particle]);
	}

//	for( int i = 0; i < 1500; i++){ phi[i] = 0 ;}
	// random particle getnerated.

	double psi_1=0, cosphi=0;
	double psi_scan=0;

// 1st my idea psi = phi's weighted average

	for(int i=0; i<N_particle; i++){
	cosphi += TMath::Cos(phi[i]);
	}
	psi_1 = TMath::ACos(cosphi/N_particle); 
	cout << psi_1 << endl;
	TLine *l = new TLine(psi_1,0, psi_1,1.5);
	l->SetLineColor(kViolet);
	l->SetLineWidth(2);
	c1->cd(1);
	h1->Draw();
	l->Draw();	
// 2nd 0~pi scan
	c1->cd(2);
	h2->Draw();
	int j=0;	
	double scan_max=0;
	double scan_x=0;
	for(psi_scan=0; psi_scan<Pi; psi_scan = psi_scan+ 1./100*Pi){
		cosphi=0;
		j=j+1;
		for (i =0; i<N_particle; i++){
			double dphi = phi[i] - psi_scan;
			cosphi += TMath::Cos(dphi);
			}	
		
		g1->SetPoint(j, psi_scan, cosphi/N_particle);
		if (scan_max < cosphi/N_particle) {scan_max = cosphi/N_particle; 
			scan_x = psi_scan;	
			}

		}
	g1->SetMarkerStyle(22);
	g1->SetMarkerSize(0.5);
	g1->SetMarkerColor(kBlue);
	double psi_2 = scan_x;
	cout << psi_2 << endl;
	TLine *l2 =  new TLine(psi_2, 0, psi_2, 1.5);
	l2->SetLineColor(kGreen);
	g1->Draw("same, p");
	l->Draw("same");
	l2->Draw("same");
}
